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Abstract 

This paper is divided in two parts. In the first part, a brief review of a spectral element 
method for the numerical solution of the incompressible Navier-Stokes equations is 
given. The method is then extended to compute buoyant flows described by the Boussi- 
nesq approximation. Free convection in closed two-dimensional cavities are computed 
and the results are in very good agreement with the available reference solutions. 
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^ ; INTRODUCTION 

At the previous conference (Wasberg et ah, 20TJT] ) we presented a method for the nu- 
merical solution of the incompressible Navier-Stokes equations by spectral element 
methods. In this contribution we will present the current status of this development 
, Qh | effort, in particular new developments to enable computation of buoyancy-driven flows 

described by the Boussinesq approximation. We will first give a review of the sta- 
tus of the basic Navier-Stokes solver. Then we describe the Boussinesq model and 
discuss its numerical implementation. Simulations of free-convection flows in closed 
two-dimensional cavities show very good agreement with available reference solutions. 

in 
o 

£2 SPECTRAL ELEMENT METHOD 



> 



O 

o 



The Navier-Stokes equations for a constant-density incompressible fluid are 



^ + u-Vu = -S7p + uV 2 u, intl, (la) 

! V • u = 0, in Q, 

with the following initial- and boundary conditions: 

u(x, 0) = u (x), x G O, V • u = 

u(x, t) — u„(x, t), x G dQ v , (lb) 
Vu(x,t) • n = 0, x e dQ . 

The boundary <9fi is divided into two parts d£l v and <9fi D , with Dirichlet velocity in- 
flow and homogeneous Neumann outflow conditions, respectively, and n is the outward 
pointing normal vector at the boundary. 

In dU), O G R d is the computational domain in d spatial dimensions, u = u(x, t) 
is the d-dimensional velocity vector, p = p(x, t) is the kinematic pressure, and v is the 
kinematic viscosity coefficient. 

To solve © we employ an implicit-explicit time splitting in which we represent the 
advective term explicitly, while we treat the diffusive term, the pressure term, and the 
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divergence equation implicitly. After semi-discretisation in time we can write ([T]) in the 
form 



(al - vV 2 )u n+1 = Vp + /(u n , u"- 1 , . . . ), (2a) 
V • u n+1 = 0, (2b) 

in which the explicit treatment of the advection term is included in the source term /. 
In the actual implementation we use the BDF2 formula for the transient term, 

TP = + O (At 2 , 

dt 2At K h 

which gives a = 3/2 At in ©, while we compute the advective contributions according 



to the operator- integration-factor (OIF) method (Maday et al, 1990 ). 



The spatial discretisation is based on a spectral element method dPatera, 1984| ); the 
computational domain is sub-divided into non-overlapping quadrilateral (in 2D) or hex- 
ahedral (in 3D) cells or elements. Within each element, a weak representation of © is 
discretised by a Galerkin method in which we choose the test and trial functions from 
bases of Legendre polynomial spaces 



G P N (x) ® P N (y) ® P N (z), (3a) 
e P N -2 0) ® P N - 2 (y) ® Pn~2 (z) . (3b) 



Note that we employ a lower order basis for the pressure spaces to avoid spurious pres- 
sure modes in the solution. The velocity variables are C 1 -continuous across element 
boundaries and are defined in the Legendre-Gauss-Lobatto points for the numerical 
integration, whereas the pressure variable is piecewise discontinuous across element 
boundaries and are defined in the interior Legendre-Gauss points. 

For the spatial discretization we now introduce the discrete Helmholtz operator, 

H = —B + uA, 
2At ' 

where A and B are the stiffness- and mass matrices in d spatial dimensions, the discrete 
divergence operator, D, and the discrete gradient operator, G. Appropriate boundary 
conditions should be included in these discrete operators. This gives the discrete equa- 
tions 

Hu n+l - Gp n+1 = Bf, (4a) 
-Du n+1 = 0, (4b) 

where the change of sign in the pressure gradient term is caused by an integration by 
parts in the construction of the weak form of the problem. This discrete system is solved 
efficiently by a second order accurate pressure correction method that can be written 

Hu* =Bf + Gp n + r, (5a) 
DQG(p n+1 - p n ) = -Du* (5b) 
u n+1 = u* + QG(p n+1 -p n ), (5c) 



where u* is an auxiliary velocity field that does not satisfy the continuity equation (l4Tt . 
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Fig. 1: Parallel speed-up for simulation of a three-dimensional fully developed 
turbulent channel flow. 

The discrete Helmholtz operator is symmetric and diagonally dominant, since the 
mass matrix of the Legendre discretisation is diagonal, and can be efficiently solved 
by the conjugate gradient method with a diagonal (Jacobi) preconditioner. The pres- 
sure operator DQG is easily computed; it is also symmetric, but ill-conditioned. The 
pressure system is solved by the preconditioned conjugate gradient method, with a mul- 
tilevel overlapping Schwarz preconditioner (Fis cher et al, 20 00). 

Earlier (Wasberg et ah, 200T| ) we presented a validation of the method for two- 
dimensional examples. Since then we have extended the method to compute the full 
three-dimensional Navier-Stokes equations. At present, turbulence simulations of a 
fully developed channel flow at Re r = 180 is in progress. Results of these simulations 
will be presented elsewhere. 

The method has good data locality and can be efficiently run on parallel computers. 
We have parallelized the code by message passing (MPI) which enables execution on 
both distributed-memory cluster and shared-memory architectures. To demonstrate the 
parallel performance, we show in Fig. [TJ the speed-up factors for a Direct Numerical 
Simulation of three-dimensional fully developed turbulent channel flow using approx- 
imately 250000 grid points. The computations were performed on a 16 processor SGI 
system. 



SOLUTION METHOD FOR BUOYANT FLOW 

The equations describing the dynamics of incompressible, viscous, buoyant flows under 



the Boussinesq approximation are 

V-u = 0, (6a) 
^ + u-Vu = - Vp + z/V 2 u + (3 (T - T rcf ) g, (6b) 
dT 

— + u-VT = aV 2 T, (6c) 



where T represents the temperature, a the thermal diffusivity, and (3 the coefficient of 
thermal expansion. The Boussinesq approximation is valid provided that the density 
variations, p(T), are small; in practice this means that that only small temperature 
deviations from the mean temperature are admitted. 

The relevant non-dimensional groups to characterize the flow are: 

• The Prandtl number Pr = v/a, 

• the Reynolds number Re = UL/u, and 

• the Rayleigh number Ra = gfiATL 3 / vol. 

Note that the Reynolds number is only relevant for problems with an imposed velocity 
scale. The free convection problems we consider below are completely determined by 
the Prandtl and Rayleigh numbers. 

In this section we will describe the solution method for the Boussinesq problem. 
Note that the buoyancy effect is accounted for in © through the solution of an addi- 
tional scalar advection-diffusion equation and an extra source term in the momentum 
equations. 

The key to efficient and accurate solution of the Boussinesq/Navier-Stokes system 
is to use an implicit-explicit splitting of diffusive and advective terms. In particular, 
if the advection/diffusion equations are solved by an implicit/explicit procedure, the 
temperature equation can be decoupled from the remaining Navier-Stokes equations, 
and the buoyancy source term can be calculated first and fed directly to the Navier- 
Stokes solver. For illustrative purposes we will discuss the solution procedure in term 
of an implicit-explicit first order Euler time discretization. Note however that higher 
order methods and operator splitting are used in the actual implementation as discussed 
for the basic Navier-Stokes solver above. A first order semi-discrete solution of the 
Boussinesq system can be written 

( K V 2 T) n+1 , (7a) 

-Vp n+1 + (/tV 2 u) n+1 + a (T n+1 - T ref ) g,(7b) 

0, (7c) 

where we have changed the ordering of the equations to emphasize that the temperature 
at the new time level, T n+1 , can be obtained from old velocity data since the advection 
term is treated explicitly. A possible solution algorithm is then self-evident: 

1 . Solve the advection-diffusion equation to obtain T n+1 . 

2. Calculate the buoyancy source term. 
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3. Calculate the explicit contributions to the momentum equations. 



Fig. 2: Streamlines and temperature distribution for the side-heated buoyant cavity 
flow at Ra = 10 3 . 



4. Solve the remaining Stokes problem with the pressure correction method to ob- 
tain u n+1 andp n+1 . 

The actual implementation is based on higher-order methods and the operator in- 
tegration factor splitting method described above for the advection/diffusion equations 
(both for temperature and momentum). The method uses second order accurate inte- 
grators; for the advection terms we use an adaptive Runge-Kutta method, while the 
implicit parts are solved by the second order implicit Euler scheme (BDF2). 

SIMULATIONS OF FREE CONVECTION CAVITIES 

We have performed simulation of the free convection in two-dimensional square and 
rectangular cavities. Cavity flows are often used as test cases for code validation, be- 
cause they are simple to set up and reliable reference solutions are readily available. 
Furthermore, thermal cavity flows display a plethora of interesting fluid dynamic phe- 
nomena, and they are important prototype flows for a wide range of practical techno- 
logical problems, including ventilation, crystal growth in liquids, nuclear reactor safety, 
and the design of high-powered laser systems. 

Differentially heated square cavity 

The steady-state differentially heated square cavity flow was the subject of one of the 



first benchmark comparison exercises reported in (de Vahl Davis & Jones, 1983). The 
benchmark results produced in that exercise are given in (de Vahl Davis, 1983| ). The 
results of de Vahl Davis were produced, for Rayleigh numbers in the range 10 3 -10 6 , 
using a stream-function/vorticity formulation discretised by a second-order finite dif- 
ference method on a regular mesh. Later, more accurate results obtained by a sec- 
ond order finite volume method on higher resolution non-uniform grids were presented 



in (Hortmann et ah, 1990). 



The problem comprises a square box of side length L x = L y = L filled with a 
Boussinesq fluid characterized by a Prandtl number, Pr = 0.71. The vertical walls are 
kept at a constant temperature, T hot and T co i d , respectively while the horizontal lid and 
bottom are insulated with zero heat flux. The direction of gravity is downwards, i.e. in 
the negative y-direction. 

We performed calculations for 10 3 < Ra < 10 6 , and at these Rayleigh numbers the 
flow is stationary. We show the computed flow field and temperature distributions in 
Figs.EHU 




Fig. 3: Streamlines and temperature distribution for the side-heated buoyant cavity 
flow at Ra = 10 4 . 




Fig. 4: Streamlines and temperature distribution for the side-heated buoyant cavity 
flow at Ra = 10 5 . 




Fig. 5: Streamlines and temperature distribution for the side-heated buoyant cavity 
flow at Ra = 10 6 . 



Table 1: Computed Nusselt numbers for the square cavity compared to the extrapolated 
results of the reference solutions of de Vahl Davis (1983) and Hortmann et al. (1990) 



Rayleighno: 10 4 10 5 10 6 

Present results: 2.245 4.522 8.825 

de Vahl Davis: 2.243 4.519 8.800 

Hortmann et al: 2.245 4.522 8.825 



The most important diagnostic connected to the free convection cavity flow is the 
average Nusselt number, which expresses the non-dimensional heat flux across the cav- 
ity. The Nusselt number is usually calculated at a vertical line, typically the hot and 
the cold wall. For consistency with the weak Galerkin formulation, we have however 
chosen to compute a global Nusselt number given by 

Nu = (8a) 

where Q is the calculated global heat flux through the cavity 

Q = / uT — a--—dxdy, (8b) 
Jo Jo dx 

and the reference value, Q , is the corresponding heat flux if the heat transfer were by 
pure conduction 

Q = L x L y ^—- = L y aAT. (8c) 

We have confirmed that the computed values of the average global Nusselt number does 
indeed agree with the average wall Nusselt numbers. 

We performed simulations using M = 4 x 4 elements varying the resolution in each 
element from N = 6 x 6 to N = 24 x 24. In Figs.|6]-|8]we show the grid convergence 
of the computed Nusselt numbers compared to the previously reported benchmark re- 
sults (de Vah l Davis, 1983| ) and ([Hortma nn et al, 1990| ). Note the excellent agreement 
with the reference data; even the coarsest resolution (i.e. 24 x 24) produces solutions 
that are essentially converged except at the highest Rayleigh number. In Table [l] we 
compare the Nusselt numbers obtained at the finest grid with the 'grid-independent' 
values from the reference solutions obtained by Richardson extrapolation. 

Bottom-heated square cavity 

It is also interesting to consider the case in which the cavity is heated from below 
instead of from the vertical walls as in the above examples. When the heated walls are 
aligned with the direction of gravity the circulation, and hence convection, is set up at 
small Rayleigh numbers. Although the bottom-heated case corresponds to a genuinely 
unstable situation; gravity will act against the instability caused by the temperature 
difference and produce a regime of pure conduction at low Ra. To illustrate this we 
show the effect of conduction expressed by the deviatory Nusselt number, i.e. Nu — 
1, for the two cases in Fig. |5| Note that whereas there is a smooth transition into 
the convection regime in the wall-heated case, the bottom-heated case shows a sharp 
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Fig. 6: Grid convergence of the average Nusselt number for the differentially heated 
buoyant cavity flow at Ra = 10 4 . 
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Fig. 7: Grid convergence of the average Nusselt number for the differentially heated 
buoyant cavity flow at Ra = 10 5 . 
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Fig. 8: Grid convergence of the average Nusselt number for the differentially heated 
buoyant cavity flow at Ra = 10 6 . 



transition point below which heat transfer is purely by conduction. Above the critical 
point the difference between the two cases, both with respect to the heat transfer and to 
the flow field, is small as we can see in Fig. [TO! 

Simulation of a tall cavity 

Christon et al. (120021) summarises the results of a workshop discussing the free con- 
vection in a tall cavity with aspect ratio 8:1. The comparison was performed for 
a Rayleigh number Ra = 3.4 x 10 5 , which is slightly above the transition point 
from steady-state to time-dependent flow at Ra ~ 3.1 x 10 5 . A total of 31 solutions 
were submitted to the workshop, of these a pseudo-spectral solution using 48 x 180 
modes ( Xin & Le Quere, 2002| ) was selected as the reference solution. 

We have computed the solution to this case with roughly the same resolution as 
in the steady-state computations of the square cavity reported above. We show the 
time history of the global Nusselt number in Fig.[TT] and note that the flow reaches a 
statistically steady state after approximately 1500 non-dimensional time units 
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Note that there appears to be good agreement with the reference solution mean. This is 
confirmed in Table[2]in which we give time averages of the computed Nusselt number, 
the velocity metric 
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Fig. 9: Deviatory Nusselt number for the square cavity heated from below. 




Fig. 10: Streamlines and temperature distribution for the bottom-heated buoyant 
cavity flow at Ra = 10 4 . 



Table 2: Computed average Nusselt numbers, average velocity norm, average vorticity 
norm, and oscillation period for the tall cavity compared to the reference solutions. 



M 



N 



At 



Nu 



U 



r 



4 x 20 6 x 6 6.92 x 10" 3 4.58356 0.2420 3.0332 3.404 
4 x 20 10 x 10 6.92 x 10~ 3 4.57951 0.2395 3.0172 3.411 
4 x 20 14 x 14 6.92 x 10" 4 4.58396 0.2421 3.0345 3.403 
4 x 20 14 x 14 1.38 x 10~ 3 4.58393 0.2421 3.0344 3.397 
4 x 20 14 x 14 2.76 x 10~ 3 4.58397 0.2421 3.0342 3.403 
Reference solutions: 4.57946 0.2397f 2.9998| 3.412 

f: The average velocity and vorticity were not given in (Xin & Le Quere, 2002) 
the reference values are the average of 29 solutions presented at the work- 
shop dChriston et al, 2002 1 



the vorticity metric 



and the oscillation period compared to the reference solutions. 
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Fig. 11: Time history for the Nusselt number in the differentially heated tall cavity at 

Ra = 3.4 x 10 5 . 
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